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ABSTRACT 

Chebyshev pseudospectral methods are used to compute two-dimensional 
smooth compressible flows* Grid refinement tests show that spectral accuracy 
can be obtained. Filtering is not needed if resolution is sufficiently high 
and if boundary conditions are carefully prescribed. 
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1 . Introduction 


Pseudospectral approximations to the compressible Euler equations have 
recently been studied as an alternative to second or fourth order finite 
differences [1], [2]. The motivation is to obtain the superior accuracy 
characteristic of spectral solutions to smooth incompressible flows. For 
simple linear hyperbolic models it is easy to demonstrate that spectral 
approximations are indeed far superior to finite difference approximations. 
(See, for example, Turkel [1], Hussalnl, Salas and Zang [3].) 

Hussaini, et. al. [4] have shown that spectral approximations do work in 
a variety of compressible flows. However, it has not been demonstrated that 
spectral accuracy is obtained in actual problems. Typically, more complicated 
flows have shocks and shock capturing spectral approximations introduce global 
oscillations which must be smoothed. In such cases a definite advantage in 
accuracy or convergence rate of spectral over finite difference approximations 
has not yet been established [4], [5]. For this reason we examine the use of 
spectral methods in conjunction with shock fitting algorithms. In the region 
where the solution is smooth there is hope of obtaining spectral accuracy. 

In this paper we examine the accuracy of spectral methods when applied to 
smooth compressible flows. Even in such cases, spectral solutions sometimes 
exhibit ’wiggles'*. We look at the need for smoothing the solutions to such 
problems. Two benchmark problems which are non-trivial and two-dimensional 
are used. The first is the interaction of a plane wave with a shock. The 
second is the classical Rlngleb flow. 
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2. Pseudospectral (tethod 

The novelty of the pseudospectral methods is that the solution defined at 
each grid point is represented by a global high order interpolating 
polynomial. Derivatives, when computed from the interpolant, couple all the 
points. For periodic problems Fourier interpolation is appropriate. Boundary 
value problems can use Chebyshev polynominals . Computation of the 
coefficients is done efficiently through the use of the Fast Fourier 
Transform. For continuous solutions the interpolation error decays faster 
than any polynominal of the number of mesh points. 

The problems in the next two sections use the Euler equations in non- 
conservation form 

Qt + BQx + CQy = 0 X, Ye [0,1], t > 0 (1) 

where Q is the column vector of the unknowns and B(Q), C(Q) are square 
matrices. For the shock/plane wave interaction problem of Section 3, the Y 
direction is periodic. Q is approximated by a Chebyshev— Fourier expansion 

M N/2-1 . 

Q(X,Y,t) = I I Q (tn (2) 

p=0 q=-N/2 P'* P 

where ? = 2X-1. The coefficients Qp^ are products of the Chebyshev and 
Fourier coefficients computed from the values of Q at the mesh points. The 
derivatives of the interpolant are 

^X = 2 I I Q^q’°^(t)T 
p=0 q=-N/2 P'^ P 

^2TTiqY 


M N/2-1 .. 

^Y = 2^ I I 

p=0 q=-N/2 P^ 


(t)Tp(O 


(4) 
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The coefficients are computed with the standard recursion formula [3] 

pq 

and 


p(O.l) 

pq 



(5) 


For the Ringleb problem a double Chebyshev approximation is used and the 
solution is approximated by 

M N 

Q(X,Y,t) = I I Qnn^t)T (OT (n). (6) 

p=0 q=0 ^ ^ 

Where n = 2Y-1. The derivatives in both directions are evaluated in a manner 
analogous to eq. (3). 

While the approximation of derivatives at boundaries often requires 
points outside the mesh, this is not the case for the spectral approximations. 
The derivatives use only points within the mesh and hence do not require 
special treatment . 

The time discretization used is the second order modified Euler. Let 
L(Q) denote the spatial discretization of t = nAt. Then 


Q = [l-AtL"]Q" 


(7) 


q"+ 1 .l/2[q"+(l-ltl)q] 


where 


L = L(Q). 


3. Shock/Plane Wave Interaction 

The first benchmark problem is the time-dependent interaction of a plane 
wave with an infinite normal shock. We use this to demonstrate the appearance 
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of wiggles in a case where the relevant features are not resolved. A detailed 
discussion of the problem and a comparision of finite difference computations 
with linear theory predictions can be found in Zang, et. al. [6]. We comment 
here only that for low amplitude waves whose wave fronts are nearly parallel 
to the moving shock the linear theoretical solutions are quite accurate. 

Let Xg(y,t) denote the position of an infinite shock moving from left 
to right into a gas which is quiescent except for a specified pressure wave of 
amplitude, A(x) . We allow the amplitude to vary smoothly from zero to a 

constant value so that the shock interacts with a smooth perturbation. In the 

absence of the pressure wave the shock would remain plane and move with a 

shock Mach number 

s 

The computational domain lies between some arbitrarily chosen left 
boundary x^ and the shock on the right. The y direction is 
periodic, — ^ y ^ This domain is mapped to the unit square by 

X = (x-Xj^)/(xg-Xj^) , Y = yly^ (8) 

T 

where y^ is the period in y. The dependent variables are 0 = (P u v S) 

where P is the logarithm of the pressure, u and v are the velocities in 
X and y, and S is the entropy divided by the specific heat at constant 

volume. 

The boundary conditions at Y = 0 and 1 are periodic. The right side 
is hounded by the moving shock and a shock fitting algorithm is used to 
determine the flow variables and move the shock. The left boundary is 

supersonic Inflow so all variables are specified. 

Table I shows the RMS error for the acoustic transmission coefficient of 
an incident 10^ pressure wave for A = 0.001 with an Mg = 3 shock on three 
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2 W 2 

different Chebyshev grids. The error is defined by e = (A^-A^) /M. The 

p=l P 

transmission coefficient k' is taken as the fundamental Fourier amplitude 

P 

computed by a Fourier transform in the Y direction at each grid point in 

X. The linearized solution is k ' . 

e 


TABLE I 

Rms Error in Acoustic Transmission Coefficient 


Number of Chebyshev Modes 

RMS Error 

8 

13.0 

16 

2.4 

32 

0.062 


Figure 1 shows k' as a function of x behind the shock for the N = 
16 and 32 Chebyshev grids. The solid line shows the linear theory results 
for a constant amplitude wave. The numerical wave is started up smoothly. 
Because the 16 point mesh cannot resolve the startup rise, large oscillations 
are present. (Note that at the time chosen the beginning of the wave occurs 
near the coarsest grid spacing.) If the solution is adequately resolved as in 
the 32 point calculation, the oscillations are almost eliminated. 

For this problem it has not been necessary to smooth or filter for the 
purpose of stability. Aesthetic reasons may lead one to cosmetically filter 
the solution. However, it must be remembered that these oscillations Indicate 
that the solution is not adequately resolved. See also Gottlieb et. al. [7] 
who show that stability can be obtained without filtering if resolution is 


adequate. 
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X 

Figure 1: Acoustic transmission coefficient computed with a 16 (circles) 

and 32 (diamonds) point Chebyshev grid. The solid line is the 
linear theory prediction. 


4. Rlngleb Flow 

The classical Ringleb flow is used for the second benchmark problem. We 
use this flow to test the algorithm on a smooth, steady, two-dimensional 
supersonic and transonic problem for which an exact solution is available. We 
also use it to bring out aspects associated with the specification of boundary 
conditions for Chebyshev methods. 

The dependent variables Q = (P u v)^ were chosen where (u,v) are the 
Cartesian velocity components. Since the flow is homentropic we simply set 
the entropy constant. The flow geometry is computed from the exact solution 
by specifying the Mach numbers at the Inflow and outflow along the lower 
streamline and the outflow Mach number along the upper streamline. Figure 2 
shows the Mach contours of the transonic problem used. The geometry is mapped 
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onto a square by a transformation to ((|>,4^), the potential-streamfunction 
coordinates, which are computed from the exact solution. A double Chebyshev 
grid is then used in these coordinates. Figure 3 shows a 17 x 9 point grid 
for the flow of fig. 2. The supersonic flow uses only the exit portion of the 
channel. 



Figure 2: Mach number contours 

for the exact solution 
to the Ringleb flow. 



Figure 3: The 17 x 9 Chebyshe- 

Chebyshev grid for the 
flow in fig. 2. 


In the mapped coordinate system correspond to (X,Y) in eq. 

(1). The coefficient matrices are 



u 

•e- 

Y<t) 

y 


V 

Y'l' 

X 

Yt 

y 

B = 

2 

£_<i 
Y X 

u 

0 

c = 

2 

Y X 

V 

0 


2 

0 

u 


2 

0 

V 


Y y 




_ Y y 




B 


( 11 ) 
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where a is the sound speed and y is the ratio of specific heats. The 
contravariant velocity components are 


U = u(b + vd) 

X 

( 12 ) 


V = + vip 


The specification of the boundary conditions has turned out to be a most 
important aspect of computing the Ringleb problem. Reference [7] details 
several studies applied to finite difference methods. We now describe the 
approach which works best with the Chebyshev methods. 

The Ringleb problem requires several types of boundary conditions. 
First, the upper and lower streamlines (ab and cd in fig. 3) are treated as 
impermeable boundaries, hereafter referred to as walls. The outflow 
boundary bd is chosen to be supersonic. Finally, the inflow boundary ac 
can be either a subsonic or supersonic boundary, depending on where it is 
placed along the channel. 

For the wall boundaries the tangential momentum equation can be written 
as 

2 

^t ^^Vx'^Vy^ ^ ° (13) 

The equation is left in this form without explicitly writing the 4> 
derivative of the contravariant velocity, U, because the derivatives of u 
and V are available from the Chebyshev interpolant. The spatial derivatives 
directly at the wall are computed as described in Section 2. The time 
discretization is performed as in eq. (7). From the fact that the 
contravariant velocity, V = 0 along the wall, u and v can be determined. 
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Particular care must be used in specifying the wall pressure when using 
spectral methods. An example of the disastrous results which can occur when 
boundary conditions are overspecified can be seen in reference [4]. Computing 
the pressure from the enthalpy or directly from the pressure equation are also 
unsatisfactory. Such boundary conditions produce wiggles even for finite 
difference computations. The only approach which works effectively is to use 
the compatibility relation for the characteristics intersecting the wall from 
the Interior of the flow. 

By combining the pressure and normal momentum equations, an equation for 
the pressure is 

= + a|vi|)|p, - fup.+rfu.i}; +u.<}) +V, 1 }^ +v.<J) ) 
t 4' '-<})''• <1) X i}' y <J> y-* 

(14) 

where the upper sign applies to the lower wall and the lower sign to the upper 
wall boundary. The spatial derivatives are again computed from the Chebyshev 
interpolant and no special treatment is needed. The equation is updated 
according to eq. (7). 

The supersonic outflow and inflow boundaries pose no difficulties. At 
the inflow all the quantities are specified. The outflow requires no boundary 
condition, either physical or numerical. Unlike typical finite difference 
methods, particularly high order ones, the Chebyshev discretization does not 
require any so-called ''numerical" boundary conditions. 

Finally, for the subsonic inflow we specify the total enthalpy and the 
angle of the flow. Typically this leads to a faster approach to the steady 
state. A compatibility condition combining the normal momentum equation and 
the pressure equation is 
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(U-a|7*|)P^ - -^[u^+n(u^+^+v^+p] 


= - y(u,i|; +u.(j) +v,\|; +v.<l) 1 
'■ <|)^x (j)^y'' 


(15) 


Since the total enthalpy is taken to be a constant along the inflow boundary, 
another relation between P and U can be obtained by differentiating the 
total enthalpy equation in time 


P 


t 


U 

|V(|)|^ 



(16) 


Solving eq. (14) and (15) allows both and to be computed. They too 

are updated according to eq. (7). From the computed U and the fact that 
V = 0, the Cartesian velocities are calculated. 


TABLE II 

Maximum Error in p for MacCormack and Spectral Computation 
of Supersonic Ringleb Flow 


Grid 

MacCormack 

Spectral 

5x5 

2.2 X 10“2 

7.5 X 10“'^ 

9x9 

4.1 X 10"3 

1.1 X 10“^ 

17 X 17 

1.0 X 10"^ 

6.6 X 10”^^ 


The fully supersonic flow is a relatively easy problem to compute. A 
9x9 grid used is shown in fig. 4. The exact solution was chosen as the 
initial condition and the computations were run long enough for errors to 
propagate out of the mesh. The time steps were kept small so that the errors 
would be dominated by the spatial discretization. For the 9x9 computation 
2000 time steps were used. The Mach contours for that solution are shown In 
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fig. 5. A grid refinement study is presented in Table II where the maximum 
error in the pressure from the spectral calculations are compared to second 
order MacCormack finite difference results. The superior error convergence 
for the spectral computations is clear. 



Chebyshev grid used 
for the supersonic 
computation. 



Figure 5: Computed Mach number 

contours of the super- 
sonic flow for the 
9x9 grid. 


The computation of the transonic flow depicted in fig. 2 is more 
difficult than the supersonic flow of fig. 5. The reason is the presence of 
the sonic line and the rapid expansion to sonic conditions along the inner 
wall. The computations were started with the exact solution and run for 
approximately the same length of physical time. The slow, explicit time 
integration method used does not allow relaxation to convergence. The Mach 
contours of a 17 x 9 point calculation are shown in fig. 6 and can be compared 
directly with fig. 2. The largest errors occur near the high curvature 
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section of the lower wall near the sonic line, and at the lower inflow 
corner. A grid refinement study is shown in Table III. Though the results 
are not as spectacular as the supersonic case, the spectral still outperforms 
the finite difference computations. 

Finally, no filtering was needed for the Ringleb problem either for the 
supersonic or the transonic cases. Solutions with wiggles result from 
boundary conditions other than the ones which we described. Application of 
the compatibility relations at the boundaries appears to be the best approach. 


TABLE III 

Maximum Error in p for MacCormack and Spectral Computation 
of Transonic Ringleb Flow 


Grid 

MacCormack 

Spectral 

9x5 

2.6 X 10“2 

2.2 X 10~2 

17 X 9 

1.1 X 10“2 

1.9 X 10"^ 

33 X 17 

3.2 X 10"^ 

5.0 X 10“5 
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5 . CONCLUSIONS 

We have shown that for two-dimensional smooth flows it is possible to 
obtain spectral accuracy characteristic of more simple problems. The first 
problem considered, that of the shock/plane wave interaction, needed no 
smoothing for stability. Oscillations were significant only if the flow was 
not well resolved. The Ringleb problem provided a more general boundary value 
test. Careful specification of the boundary conditions allowed the 
computation to be performed without smoothing. 

Though the spectral method is superior to the finite difference method in 
terms of accuracy, the computation times for the spectral are far longer. One 
major difficulty comes from the use of an explicit time differencing 
procedure. Even for finite difference computations convergence to a steady 
state is very slow without some acceleration procedure. The spectral 
computations have the added disadvantage that the time step, which depends on 
the grid spacing, varies with rather than N. The widespread use of the 
spectral methods will even more strongly depend on the development of fast 
relaxation methods for the Euler equations. 
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